#####################################################
### Code to create Figure 6, Appendix Tables 9-11 ###
#####################################################



################################################################################
### Setup 

# Clear workspace, call packages
rm(list=ls())
library(MASS)
library(ggplot2)
library(dplyr)
library(stargazer)
library(sjPlot)
library(jtools)
library(lfe)
library(data.table)
library(lubridate)
library(Hmisc)
library(abmisc)
library(tidyverse)


### Merge member covariates

# Load NOMINATE data
nominate <- read.csv("HSall_members2021.csv")

# Load member lists with GovTrack data
c111 <- read.csv("Members by Congress 111th.csv")
c112 <- read.csv("Members by Congress 112th.csv")
c113 <- read.csv("Members by Congress 113th.csv")
c114 <- read.csv("Members by Congress 114th.csv")
c115 <- read.csv("Members by Congress 115th.csv")
c116 <- read.csv("Members by Congress 116th.csv")

# Push GovTrack to numeric
c111$GovTrack <- as.numeric(c111$GovTrack)
c112$GovTrack <- as.numeric(c112$GovTrack)
c113$GovTrack <- as.numeric(c113$GovTrack)
c114$GovTrack <- as.numeric(c114$GovTrack)
c115$GovTrack <- as.numeric(c115$GovTrack)
c116$GovTrack <- as.numeric(c116$GovTrack)

# Create distance metric from party center
c111$govdist <- abs(c111$GovTrack - mean(c111$GovTrack,na.rm=TRUE))
c112$govdist <- abs(c112$GovTrack - mean(c112$GovTrack,na.rm=TRUE))
c113$govdist <- abs(c113$GovTrack - mean(c113$GovTrack,na.rm=TRUE))
c114$govdist <- abs(c114$GovTrack - mean(c114$GovTrack,na.rm=TRUE))
c115$govdist <- abs(c115$GovTrack - mean(c115$GovTrack,na.rm=TRUE))
c116$govdist <- abs(c116$GovTrack - mean(c116$GovTrack,na.rm=TRUE))

# Split NOMINATE data into congresses to merge with GovTrack data
nominate111 <- dplyr::select(filter(nominate, congress == "111"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)
nominate112 <- dplyr::select(filter(nominate, congress == "112"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)
nominate113 <- dplyr::select(filter(nominate, congress == "113"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)
nominate114 <- dplyr::select(filter(nominate, congress == "114"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)
nominate115 <- dplyr::select(filter(nominate, congress == "115"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)
nominate116 <- dplyr::select(filter(nominate, congress == "116"), 
                             icpsr, chamber, nominate_dim1, nokken_poole_dim1)

# Merge GovTrack and NOMINATE data
m111 <- merge(c111, nominate111, by = "icpsr")
m112 <- merge(c112, nominate112, by = "icpsr")
m113 <- merge(c113, nominate113, by = "icpsr")
m114 <- merge(c114, nominate114, by = "icpsr")
m115 <- merge(c115, nominate115, by = "icpsr")
m116 <- merge(c116, nominate116, by = "icpsr")

# Create NOMINATE distance measure
m111$DWdist <- abs(m111$nominate_dim1 - mean(m111$nominate_dim1))
m112$DWdist <- abs(m112$nominate_dim1 - mean(m112$nominate_dim1))
m113$DWdist <- abs(m113$nominate_dim1 - mean(m113$nominate_dim1))
m114$DWdist <- abs(m114$nominate_dim1 - mean(m114$nominate_dim1))
m115$DWdist <- abs(m115$nominate_dim1 - mean(m115$nominate_dim1))
m116$DWdist <- abs(m116$nominate_dim1 - mean(m116$nominate_dim1))

# Combine congresses
combined <- rbind(m111, m112, m113, m114, m115, m116)

# Create year-name variable for merging
combined$icpsrcongress <- paste0(as.character(combined$icpsr),"-", as.character(combined$congress))


### Format tweets data

# Load tweets, filter out RTs
tweets <- fread("activetweets111-116.csv") %>% 
  filter(RT==0)

# Get number of polarizing tweets and number of tweets by member-congress
congresstweets <- tweets %>%
  group_by(twitter_lower, congress) %>%
  dplyr::summarize(
    tweets = n(),
    polarizing = sum(polarizing)
  ) %>%
  ungroup()

# Load handles
handles <- read.csv("US Congress Handles Master List.csv")
handles <- dplyr::select(handles, icpsr, twitter_lower, Official)

# Merge tweets and handles (adds ICPSR to tweets)
congresstweets <- merge(congresstweets, handles, by = "twitter_lower")

# Create percent polarizing
congresstweets$pct.polarizing <- ((congresstweets$polarizing / congresstweets$tweets))

# Create icpsrcongress to match with combined covariate data
congresstweets$icpsrcongress <-  paste0(as.character(congresstweets$icpsr),"-", as.character(congresstweets$congress))

# Merge tweets data with covariate data
congresstweets <- merge(congresstweets, combined, by = "icpsrcongress")

# Dummy variable for Republican
congresstweets$republican <- ifelse(congresstweets$party == "R", 1, 0)



################################################################################
### Percentage Polarizing Models


### Run models

# Base models
govmodel <- felm(pct.polarizing ~ govdist + Pres.Party + PVIABS + chamber.y + 
                   female + republican | congress.x, data = congresstweets)
DWmodel <- felm(pct.polarizing ~ DWdist + Pres.Party + PVIABS + chamber.y + 
                  female + republican| congress.x, data = congresstweets)

### Pre/post-trump split models
govmodelpre <- felm(pct.polarizing ~ govdist + Pres.Party + PVIABS + chamber.y + 
                      female | congress.x, 
                    data = filter(congresstweets, congress.x < 115))
govmodelpost <- felm(pct.polarizing ~ govdist + Pres.Party + PVIABS + chamber.y + 
                       female | congress.x, 
                     data = filter(congresstweets, congress.x > 114))

DWmodelpre <- felm(pct.polarizing ~ DWdist + Pres.Party + PVIABS + chamber.y + 
                     female | congress.x, 
                   data = filter(congresstweets, congress.x < 115))
DWmodelpost <- felm(pct.polarizing ~ DWdist + Pres.Party + PVIABS + chamber.y + 
                      female | congress.x, 
                    data = filter(congresstweets, congress.x > 114))


### Model Appendix Tables 9 and 11
stargazer(govmodel, DWmodel, type="latex",
          covariate.labels = c("Ideological Extremity (GovTrack)", 
                               "Ideological Extremity (DW-Nominate)", 
                               "President's Party", "|PVI|", "Senate", 
                               "Female", "Republican"),
          dep.var.labels   = ("Percent Polarizing"))


stargazer(govmodelpre, DWmodelpre, govmodelpost, DWmodelpost, type="latex", 
          covariate.labels = c("Ideological Extremity (GovTrack)", 
                               "Ideological Extremity (DW-Nominate)",
                               "President's Party", "|PVI|", "Senate", "Female"),
          dep.var.labels   = c("Percent Polarizing"),
          column.labels = c("Pre", "Post"))



#################################################################################
### Plot results (Figure 6)

### Data for effect sizes

# Flexible function for each plot
effPlotDat <- function(models, plot_vars, mod_df=congresstweets){
  if(length(models)==1){
    mod <- models[[1]]
    # Pull coefficients
    coef_tab <- summary(mod)$coefficients
    coefs <- coef_tab[plot_vars,'Estimate']
    cis <- coef_tab[plot_vars,'Std. Error']*1.96
    
    # Put in terms of 2 SD increase for continuous variables
    for(ii in seq_along(plot_vars)){
      if(plot_vars[ii] != 'Pres.Party'){
        coefs[ii] <- coefs[ii]*sd(mod_df[,plot_vars[ii]], na.rm=TRUE)*2
        cis[ii] <- cis[ii]*sd(mod_df[,plot_vars[ii]], na.rm=TRUE)*2
      }
    }
    
    # Combine/return
    effs <- data.frame(coef=coefs, ci=cis)
    effs$coef[2] <- abs(effs$coef[2]) # because the President's party variable is coded as being *in* the president's party
    
  } else {
    effs <- sapply(models, function(mod){
      # Pull coefficients
      iter_vars <- intersect(plot_vars, rownames(mod$coefficients))
      coef_tab <- summary(mod)$coefficients
      coefs <- coef_tab[iter_vars,'Estimate']
      cis <- coef_tab[iter_vars,'Std. Error']*1.96
      
      # Put in terms of 2 SD increase for continuous variables
      for(ii in seq_along(iter_vars)){
        if(iter_vars[ii] != 'Pres.Party'){
          coefs[ii] <- coefs[ii]*sd(mod_df[,iter_vars[ii]], na.rm=TRUE)*2
          cis[ii] <- cis[ii]*sd(mod_df[,iter_vars[ii]], na.rm=TRUE)*2
        }
      }
      
      # Combine/return
      eff <- data.frame(coef=coefs, ci=cis)
      eff
    }) %>% 
      t() %>%
      as.data.frame() %>%
      mutate(coef=num(coef), ci=num(ci))
  }
  # Deal with president's party reverse coding
  if(any(str_detect(plot_vars, 'Party'))){
    effs$coef <- abs(effs$coef)
  }
  effs
}

# President's party and electoral safety
effs_ps <- effPlotDat(models=list(govmodel), plot_vars=c('Pres.Party', 'PVIABS'))

# Two versions of ideological extremity
effs_ie <- effPlotDat(models=list(govmodel, DWmodel), plot_vars=c('govdist', 'DWdist'))

# Pre and post trump
effs_tr <- effPlotDat(models=list(DWmodelpre, DWmodelpost), plot_vars='Pres.Party')

# Plot president's party
barCenters <- barplot(c(effs_ps[1,1], effs_tr[,1]), ylim=c(0, 0.15))
jpeg('incivilityMods_effects_party.jpg', width=6, height=6, units='in', res=300)
par(mar=c(3.1, 4.1, 2.1, 1.1))
barplot(c(effs_ps[1,1], effs_tr[,1]), ylim=c(0, 0.15), 
        ylab='Effect on Proportion Polarizing', 
        names.arg=c("Overall", 'Obama', 'Trump'))
segments(x0=barCenters[1], 
         y0=effs_ps[1,1]-effs_ps[1,2], 
         y1=effs_ps[1,1]+effs_ps[1,2], 
         lwd=3)
arrows(x0=barCenters[1], 
       y0=effs_ps[1,1]-effs_ps[1,2], 
       y1=effs_ps[1,1]+effs_ps[1,2], 
       lwd=3, angle=90, code=3, length=0.05)  
text(x=barCenters[1], 
     y=effs_ps[1,1]/4, 
     paste0(round(effs_ps[1,1], 2), ' +/- ', round(effs_ps[1,2], 3)))
for(ii in 1:2){
  segments(x0=barCenters[ii+1], 
           y0=effs_tr[ii,1]-effs_tr[ii,2], 
           y1=effs_tr[ii,1]+effs_tr[ii,2], 
           lwd=3)
  arrows(x0=barCenters[ii+1], 
         y0=effs_tr[ii,1]-effs_tr[ii,2], 
         y1=effs_tr[ii,1]+effs_tr[ii,2], 
         lwd=3, angle=90, code=3, length=0.05)  
  text(x=barCenters[ii+1], 
       y=effs_tr[ii,1]/4, 
       paste0(round(effs_tr[ii,1], 2), ' +/- ', round(effs_tr[ii,2], 3)))
}
par(mar=c(5.1, 4.1, 4.1, 2.1))
dev.off()

# Electoral safety
barCenters <- barplot(effs_ps[2,1], ylim=c(0, 0.15))
jpeg('incivilityMods_effects_electoral_safety.jpg', 
     width=6, height=6, units='in', res=300)
par(mar=c(3.1, 4.1, 2.1, 1.1))
barplot(effs_ps[2,1], ylim=c(0, 0.15), xlim=c(-0.5, 2), 
        ylab='Effect on Proportion Polarizing', names.arg='')
segments(x0=barCenters, 
         y0=effs_ps[2,1]-effs_ps[2,2], 
         y1=effs_ps[2,1]+effs_ps[2,2], 
         lwd=3)
arrows(x0=barCenters, 
       y0=effs_ps[2,1]-effs_ps[2,2], 
       y1=effs_ps[2,1]+effs_ps[2,2], 
       lwd=3, angle=90, code=3, length=0.05)  
text(x=barCenters, 
     y=effs_ps[2,1]/4, 
     paste0(round(effs_ps[2,1], 3), ' +/- ', round(effs_ps[2,2], 3)))
par(mar=c(5.1, 4.1, 4.1, 2.1))
dev.off()

# Plot two ways to measure ideological extremity
barCenters <- barplot(effs_ie[,1], ylim=c(0, 0.15))
jpeg('incivilityMods_effects_ideology.jpg', width=6, height=6, units='in', res=300)
par(mar=c(3.1, 4.1, 2.1, 1.1))
barplot(effs_ie[,1], ylim=c(0, 0.15), ylab='Effect on Proportion Polarizing', 
        names.arg=c("GovTrack", 'DW-NOMINATE'))
for(ii in 1:2){
  segments(x0=barCenters[ii], 
           y0=effs_ie[ii,1]-effs_ie[ii,2], 
           y1=effs_ie[ii,1]+effs_ie[ii,2], 
           lwd=3)
  arrows(x0=barCenters[ii], 
         y0=effs_ie[ii,1]-effs_ie[ii,2], 
         y1=effs_ie[ii,1]+effs_ie[ii,2], 
         lwd=3, angle=90, code=3, length=0.05)  
  text(x=barCenters[ii], 
       y=effs_ie[ii,1]/4, 
       paste0(round(effs_ie[ii,1], 2), ' +/- ', round(effs_ie[ii,2], 3)))
}
par(mar=c(5.1, 4.1, 4.1, 2.1))
dev.off()



################################################################################
### Official versus Campaign / Account level Models (Appendix Table 10)

# Get data by handle rather than member
congresstweets <- tweets %>%
  group_by(twitter_lower, congress) %>%
  dplyr::summarize(
    tweets = n(),
    polarizing = sum(polarizing)
  ) %>%
  ungroup()

# Load handles and merge in ICPSR for member fixed effects
handles <- read.csv("US Congress Handles Master List.csv")
handles <- dplyr::select(handles, icpsr, twitter_lower, Official)
handletweets <- merge(congresstweets, handles, by = "twitter_lower")

# Create percent polarizing and icpsrcongress
handletweets <- handletweets %>%
  mutate(pct.polarizing = polarizing/tweets,
         icpsrcongress = paste0(icpsr, '-', congress))

# Run models
govmodelhandles <- felm(pct.polarizing ~ govdist + Pres.Party + PVIABS + 
                          Chamber.Majority + chamber.y + female + republican + 
                          Official | congress.x, data = handletweets)
DWmodelhandles <- felm(pct.polarizing ~ DWdist + Pres.Party + PVIABS + 
                         Chamber.Majority + chamber.y + female + republican + 
                         Official | congress.x, data = handletweets)

# Create table
stargazer(govmodelhandles, DWmodelhandles, type='latex',
          covariate.labels = c("Ideological Extremity (GovTrack)", 
                               "Ideological Extremity (DW-Nominate)", 
                               "President's Party", "|PVI|", "Senate", 
                               "Female", "Republican", "Official"),
          dep.var.labels   = ("Percent Polarizing"))
